Frequency estimation

ABSTRACT

Determining a property of an electrical signal at a node instant at the current instant (INSTC); determining a previous function indicative of a property of the signal at the node at a previous instant; determining an error function in a previous frequency estimation at the node at a previous instant; determining a first current function indicative of the property of the signal at the node at INSTC in accordance with the determined property, the previous function, and the error function; receiving, from another node, a second current function indicative of a property of an electrical signal at the at least one other node at INSTC; combining the first current function and the second current function to produce a current combined function (CCF) indicative of a property of the signal at the node at INSTC; and estimating a current frequency of the signal at the node in accordance with the CCF.

FIELD OF INVENTION

The present disclosure relates to a method and apparatus for frequencyestimation. More specifically, but not exclusively, the presentdisclosure relates to a method and apparatus for estimating frequency ina three-phase power system.

BACKGROUND TO THE INVENTION

Electricity networks are undergoing wholesale changes both from thegeneration and the user (load) sides. Major challenges in this directionare envisaged to be the management of largely dynamic load levels, dueto e.g., charging a large number of plugin electric vehicles (PEVs), andthe duality between loads and supplies, for instance, when PEVs are usedin the “vehicle to grid” mode to mitigate power shortage and systemimbalances. Generation, historically aggregated into large power plantsand far from the user, is gradually moving towards being located at thedistribution level and based on renewable sources, that is,intrinsically intermittent. This will require enhanced flexibility ofthe grid and the ability to accommodate islanding and microgrids.

The idea behind smart distributed grids and microgrids is to balance asmuch as possible locally between production and consumption. However,the deployment of intermittent renewable sources will inevitably lead tofrequent imbalance between supply and demand, as exemplified by thedifficulties in maintaining system balance due to wind powervariability. Signal processing is certain to play a significant role indealing with the complexity and uncertainty associated with the smartgrid and in accurate real-time estimation of the system parameters, inboth balanced and unbalanced conditions.

Unexpected frequency variations from the nominal value can triggerabnormal power system conditions that can propagate and aggregate withinthe power grid system. Some of the major abnormal power systemconditions are discussed below:

-   -   Imbalance in the generation (G) and load (L). In the smart grid,        the system will frequently switch between the main grid (MG) and        microgrids (mG), with parts of the system completely switching        off the MG for prolonged periods of time (islanding). The system        frequency rises for G>L and decreases for G<L.    -   Single- and dual-phase faults. The system frequency is derived        from the relationship between the three-phase voltages. Faults        in one or two phases and voltage sags (sudden drop in voltage        for a short period of time) will cause an incorrect frequency        estimate and alarm spread through the system, although the        actual system frequency was correct.    -   Dual character of load-supply. The smart grid employs dynamic        loads and dual load-generator devices, such as PEVs, which can        give the energy back to the grid in the case of emergency.        Frequent switching will cause problems with reactive power,        whose drifting causes oscillations of power levels and harmonics        in frequency.    -   Harmonics. Some loads (power supplies, motors, heating elements)        have nonlinear voltage to current characteristics and introduce        harmonics, which may be slowly floating and not integer        multiplies of system frequency. They may cause resonance in the        system leading to significant increase in currents and        overheating of transformers. Switching on the shunt capacitors        for reactive power compensation also causes strong transients        and harmonics that are damaging to some equipment.    -   Transient stability issues. Faults and short circuits could        trigger instability, and actions such as shedding loads (or        generators) that are needed to mitigate the problem require        accurate frequency estimation.

Accurate and fast frequency tracking is therefore a prerequisite toenable the system to respond quickly to such problems. Approaches tofrequency estimation from a single phase in three-phase power networksresult in non-unique solutions. Consequently, robust frequencyestimators should consider all of the three-phase voltages. However, theprocessing required to determine the frequency from each of the threesignals in a three-phase system is computationally expensive.Consequently, various processing algorithms have been proposed in orderto attempt to reduce this computational complexity.

One such means for reducing computational complexity is the Clarke's αβtransform which produces a complex-valued signal from the three-phasevoltages, where system frequency is obtained from the phase of thissignal. Complex domain solutions for frequency tracking includephase-locked loops (PLLs), least squares, and demodulation methods.Recently, adaptive tracking algorithms based on the minimization of meansquare error (MSE) have become a standard, as they are naturally suitedto deal with noise, harmonics, and non-stationary environments. However,unbalanced events make it difficult to calculate phase angle, as in thiscase the complex-valued signal obtained from an unbalanced three-phasevoltage source is represented as an orthogonal sum of positive(reflecting the energy transfer between generators and consumers) andnegative (indicating imbalance between three-phase voltages) sequences.Standard complex linear adaptive filters can only cater for the positivesequences, whereas the negative sequences introduce a modelling errorthat oscillates at twice the system frequency. Hence, current systemsfor frequency tracking in three-phase systems are nowhere near optimal.

SUMMARY OF INVENTION

Embodiments of the present invention attempt to mitigate at least someof the above-mentioned problems.

According to an aspect of the invention a method is provided forestimating a frequency of an electrical signal at a node in anelectrical system having a plurality of nodes. The method comprisesdetermining, at a node of an electrical system having a plurality ofnodes, at least one property of an electrical signal at the node at acurrent instant in time. In addition, the method comprises determining aprevious function indicative of one or more properties of the electricalsignal at the node at a previous instant in time. The method alsocomprises determining a function of an error in a previous frequencyestimation at the node at a previous instant in time. Furthermore, themethod comprises determining a current function indicative of the one ormore properties of the electrical signal at the node at the currentinstant in time in accordance with the at least one determined propertyof the electrical signal at the node at the current instant in time, theprevious function indicative of the one or more properties of theelectrical signal at the node at a previous instant in time, and thefunction of the error in the previous frequency estimation. The methodalso comprises receiving, from at least one other node of the pluralityof nodes, a current function indicative of one or more properties of anelectrical signal at the at least one other node at the current instantin time. Then, the method comprises combining the current functionindicative of one or more properties of the electrical signal at thenode and the current function indicative of one or more properties ofthe electrical signal at the at least one other node to produce acurrent combined function indicative of one or more properties of theelectrical signal at the node at the current instant in time.Furthermore, the method comprises estimating a current frequency of theelectrical signal at the node in accordance with the current combinedfunction.

The one or more properties of the electrical signal that the previousand current functions are indicative of may include a frequency of theelectrical signal. Also, the one or more properties of the electricalsignal that the previous and current functions are indicative of mayinclude a phase associated with the electrical signal.

The determined property of the electrical signal at the node at thecurrent instant in time may be determined by measuring the property ofthe electrical signal at the node. Alternatively, the may be determinedby estimating the property of the electrical signal.

The determined property of the electrical signal may be a voltage.Furthermore, the determined property of the electrical signal may be aphase associated with the electrical signal. The electrical signal maybe a three-phase electrical signal.

The method may further comprise converting the three-phase electricalsignal into a transform-domain. The current function indicative of oneor more properties of the electrical signal at the node at the currentinstant in time may be determined in accordance with the electricalsignal in the transform-domain.

The function of an error in the previous estimated frequency may bedetermined by multiplication of a calculated error in a previousestimated frequency by a weighting. The weighting may be indicative of alevel of confidence in the calculated error. Furthermore, the functionof an error may be a function of a phase-only error.

The combining may be a linear combination process in the original or atransform domain.

The method may further comprise applying a node weighting factor to thereceived current function after receiving the current functionindicative of one or more properties of the electrical signal at the atleast one other node and before performing the step of combining. Thenode weighting factor may be indicative of a confidence in the accuracyof the received current function.

The current estimated frequency {circumflex over (f)}_(i)(k) at the nodei, may be determined by:

${{\hat{f}}_{i}(k)} = {\frac{1}{2\pi\;\Delta\; T}{\sin^{- 1}\left( {{funct}(k)} \right)}}$wherein ΔT is a sampling interval, and funct(k) is the current combinedfunction at the current instant in time, k.

The previous function indicative of one or more properties of theelectrical signal at the node at a previous instant in time, the currentfunction indicative of one or more properties of the electrical signalat the node at the current instant in time, and the function receivedfrom the at least one other node indicative of one or more properties ofthe electrical signal at the at least one other node at the currentinstant in time, each may comprise a first function and a secondfunction, wherein the first function utilises the at least onedetermined property of the electrical signal, and the second functionutilises the complex conjugate of the at least one determined propertyof the signal. Furthermore, the step of combining may comprisecombining, at the node, the first function associated with the node withthe first function associated with the at least one other node toproduce a current combined first function. The step of combining mayalso comprise combining, at the node, the second function associatedwith the node with the second function associated with the at least oneother node to produce a current combined second function. Furthermore,the step of estimating may comprise estimating a current frequency ofthe at least one electrical characteristic at the node in accordancewith the current combined first function and the current combined secondfunction.

The first function and the second function may also utilise a functionof a calculated error. Furthermore, the first function and the secondfunction may also utilise a function of the conjugate of a complextransform of the at least one determined property of the electricalsignal at the node at the previous instant in time.

The first function, h, at the current instant in time, k, may bedetermined by:h(k)=h(k−m)+μp(e(k−m))q(v*(k−m))wherein m can take any integer value, h(k−m) is a previous firstfunction, p(e(k−m)) is a function of a calculated error, μ is aweighting factor applied to the calculated error, and q(v*(k−m)) is afunction of the conjugate of a complex transform of the at least onedetermined property of the electrical signal at the node at the previousinstant in time.

The second function, g, at the current instant in time, k, may bedetermined by:g(k)=g(k−m)+μp(e(k−m))q(v(k−m))wherein m can take any integer value, g(k−m) is a previous secondfunction, p(e(k−m)) is a function of the calculated error, μ is aweighting factor applied to the calculated error at the previous instantin time, and q(v(k−m)) is a function of a complex transform of the atleast one determined property of the electrical signal at the node at aprevious instant in time.

The function of the calculated error and the function of a complextransform of the at least one determined property of the electricalsignal at the node, may be determined in accordance with one or thefollowing algorithms:

${{{NCLMS}\mspace{14mu}{algorithm}\text{:}\mspace{14mu}{p\left( {e(k)} \right)}} = {e(k)}},{{q\left( {v(k)} \right)} = \frac{v(k)}{{{v(k)}}^{2}}},{or}$AP  algorithm:  p(e(k)) = e(k), q(V(k)) = V(k)(V^(T)(k)V^(*)(k) + δ I)⁻¹, or${{{LMP}\mspace{14mu}{algorithm}\text{:}\mspace{14mu}{p\left( {e(k)} \right)}} = {{\angle\;{v\left( {k + 1} \right)}} - {\angle\mspace{11mu}{\hat{v}\left( {k + 1} \right)}}}},{{q\left( {v(k)} \right)} = \frac{{jv}(k)}{{\hat{v}}^{*}\left( {k + 1} \right)}}$

Where {circumflex over (v)}(k+1) is the at least one property of anelectrical signal determined by estimation, e(k) is the function of anerror in the frequency estimation where defined ase(k)=v(k+1)−{circumflex over (v)}(k+1), V(k) is the input matrix,defined as:

${V(k)} = \begin{bmatrix}{{v(k)},{v\left( {k - 1} \right)},\ldots\mspace{14mu},{v\left( {k - m} \right)}} \\{{v\left( {k - 1} \right)},{v\left( {k - 2} \right)},\ldots\mspace{14mu},{v\left( {k - m - 1} \right)}} \\\vdots \\{{v\left( {k - n} \right)},{v\left( {k - n - 1} \right)},\ldots\mspace{14mu},{v\left( {k - m - n} \right)}}\end{bmatrix}$δ is a regularisation term, I is an identity matrix, ∠ is an angle of acomplex variable, and p(⋅) and q(⋅) are functions.

The step of producing the current combined first function may be definedby:

${{\underset{\_}{h}}_{i}(k)} = {\sum\limits_{l \in N_{i}}{{c\left( {i,l} \right)}{h_{l}(k)}}}$wherein h_(l)(k) is a current first function associated with a node ofthe plurality of nodes and c(i,l) is a first node weighting factorindicative of a confidence in the accuracy of a current first functionassociated with the respective node.

The step of producing the current combined second function may bedefined by:

${{\underset{\_}{g}}_{i}(k)} = {\sum\limits_{l \in N_{i}}{{c\left( {i,l} \right)}{g_{l}(k)}}}$

Wherein g_(i)(k) is a current second function associated with a node ofthe plurality of nodes and c(i,l) is a second node weighting factorindicative of a confidence in the accuracy of the current secondfunction associated with the respective node.

The current estimated frequency {circumflex over (f)}_(i)(k) at the nodemay be determined by:

${f_{i}(k)} = {\frac{1}{2{\pi\Delta}\; T}{\sin^{- 1}\left( {{Im}\left\{ {{{\underset{\_}{h}}_{i}(k)} + {{{\underset{\_}{a}}_{i}(k)}{{\underset{\_}{g}}_{i}(k)}}} \right\}} \right)}}$wherein ΔT is a sampling interval, Im{.} is the imaginary part operator,h_(i) (k) is a current combined first function, g_(i)(k) is a currentcombined second function, and

${{\underset{\_}{a}}_{i}(k)} = {\frac{{{- {j{Im}}}\left\{ {{\underset{\_}{h}}_{i}(k)} \right\}} + {j\sqrt{{{Im}^{2}\left\{ {{\underset{\_}{h}}_{i}(k)} \right\}} - {{{\underset{\_}{g}}_{i}(k)}}^{2}}}}{{\underset{\_}{g}}_{i}(k)}.}$

The current function indicative of the one or more properties of theelectrical signal at the node may also be determined in accordance witha conjugate of the at least one determined property of the electricalsignal.

The current function indicative of the one or more properties of theelectrical signal at the node at the current instant in time may bedetermined in accordance with one of the following state space models:

State Space Model 1

State equation:

$\begin{bmatrix}{h_{i}(k)} \\{g_{i}(k)} \\{h_{i}^{*}(k)} \\{g_{i}^{*}(k)}\end{bmatrix} = {\begin{bmatrix}{h_{i}\left( {k - m} \right)} \\{g_{i}\left( {k - m} \right)} \\{h_{i}^{*}\left( {k - m} \right)} \\{g_{i}^{*}\left( {k - m} \right)}\end{bmatrix} + {u_{i}(k)}}$

Observation equation for the node, i:

$\begin{bmatrix}{v_{i}(k)} \\{v_{i}^{*}(k)}\end{bmatrix} = {\begin{bmatrix}{{{v_{i}\left( {k - m} \right)}{h_{i}(k)}} + {{v_{i}^{*}\left( {k - m} \right)}{g_{i}(k)}}} \\{{{v_{i}^{*}\left( {k - m} \right)}{h_{i}^{*}(k)}} + {{v_{i}\left( {k - m} \right)}{g_{i}^{*}(k)}}}\end{bmatrix} + {n_{i}(k)}}$

State Space Model 2

State equation:

$\begin{bmatrix}{h_{i}(k)} \\{g_{i}(k)} \\{h_{i}^{*}(k)} \\{g_{i}^{*}(k)} \\{z_{i}^{*}(k)}\end{bmatrix} = {\begin{bmatrix}{h_{i}\left( {k - m} \right)} \\{g_{i}\left( {k - m} \right)} \\{{{z_{i}\left( {k - m} \right)}{h_{i}\left( {k - m} \right)}} + {{z_{i}^{*}\left( {k - m} \right)}{g_{i}\left( {k - m} \right)}}} \\{h_{i}^{*}\left( {k - m} \right)} \\{g_{i}^{*}\left( {k - m} \right)} \\{{{z_{i}^{*}\left( {k - m} \right)}{h_{i}^{*}\left( {k - m} \right)}} + {{z_{i}\left( {k - m} \right)}{g_{i}^{*}\left( {k - m} \right)}}}\end{bmatrix} + {u_{i}(k)}}$

Observation equation for a node i:

$\begin{bmatrix}{v_{i}(k)} \\{v_{i}^{*}(k)}\end{bmatrix} = {\begin{bmatrix}{z_{i}(k)} \\{z_{i}^{*}(k)}\end{bmatrix} + {n_{i}(k)}}$

State Space Model 3

State equation:

$\begin{bmatrix}{h_{i}(k)} \\{g_{i}(k)}\end{bmatrix} = {\begin{bmatrix}{h_{i}\left( {k - m} \right)} \\{g_{i}\left( {k - m} \right)}\end{bmatrix} + {u_{i}(k)}}$

Observation equation for a node i:v _(i)(k)=v _(i)(k−m)h _(i)(k)+v _(i)*(k−m)g _(i)(k)+n _(i)(k)

State Space Model 4

State equation:

$\begin{bmatrix}{h_{i}(k)} \\{g_{i}(k)} \\{z_{i}(k)}\end{bmatrix} = {\begin{bmatrix}{h_{i}\left( {k - m} \right)} \\{g_{i}\left( {k - m} \right)} \\{{{z_{i}\left( {k - m} \right)}{h_{i}\left( {k - m} \right)}} + {{z_{i}^{*}\left( {k - m} \right)}{g_{i}\left( {k - m} \right)}}}\end{bmatrix} + {u_{i}(k)}}$

Observation equation for a node i:v _(i)(k)=z _(i)(k)+n _(i)(k)

Wherein h_(i)(k−m) is the first function indicative of the one or moreproperties of the electrical signal at the node, i, at a previousinstant in time, k−m, g_(i)(k−m) is the second function indicative ofthe one or more properties of the electrical signal at the node at theprevious instant in time, k−m, z_(i)(k−m) is an estimation of the one ormore properties of the electrical signal at the previous instant intime, k−m, u_(i)(k) and n_(i)(k) are functions indicative of an error inthe previous frequency estimation, and v_(i)(k−m) is the determinedproperty of the electrical signal at the previous instant in time, k−m,at the node, i.

According to another aspect of the invention there is provided a methodfor estimating a frequency of a signal. The method comprises determininga current function indicative of one or more properties of a signal at anode at a current instant in time in accordance with: a previousfunction indicative of one or more properties of the signal at the nodeat a previous instant in time; a determined property of the signal; aconjugate of the determined property of the signal; and a function of anerror in a frequency estimation at the previous instant in time; andestimating a current frequency of the signal in accordance with thedetermined current function.

The determined property of the signal may be determined by ameasurement. Alternatively, the determined property of the signal may bedetermined by an estimation.

The signal may be an electrical signal. The determined property of thesignal may be a voltage of the signal. The signal may be a three-phaseelectrical signal. The method may comprise converting the three-phaseelectrical signal into a transform-domain. The current function may bedetermined in accordance with the electrical signal in thetransform-domain.

The function of an error in a frequency estimation at the previousinstant in time may be determined by multiplication of a function of acalculated error in the frequency estimation at the previous instant intime by a weighting. The weighting may be indicative of a level ofconfidence in the calculated error.

The function of an error may be a function of a phase-only error.

Each of the previous function and the current function indicative of oneor more properties of the signal at the node may be determined inaccordance with a first function and a second function. The firstfunction may utilise the at least one determined property of the signal,and the second function may utilise the complex conjugate of the atleast one determined property of the signal.

The current function indicative of one or more properties of the signalat the node at the current instant in time may be determined inaccordance with one of the state space equations noted above.

The method may further comprise receiving, from at least one other nodeconnected to the node, a current function indicative of one or moreproperties of a signal at the at least one other node at the currentinstant in time. The frequency may be estimated in accordance with thereceived current function indicative of one or more properties of thesignal at the at least one other node.

Prior to estimating the current frequency, the current functionindicative of one or more properties of the signal at the node and thecurrent function indicative of one or more properties of the signal atthe at least one other node may be combined in accordance with aweighting, wherein the weighting is indicative of a level of confidencein a quality of the frequency estimation at the respective node.

The current estimated frequency f_(i)(k) at the node i, may bedetermined by:

${{\hat{f}}_{i}(k)} = {\frac{1}{2{\pi\Delta}\; T}{\sin^{- 1}\left( {F(k)} \right)}}$wherein ΔT is a sampling interval, and F(k) is the current functionindicative of one or more properties of the signal at the node at acurrent instant in time.

According to another aspect of the invention an apparatus is providedfor estimating a frequency. The apparatus comprises a processing unitarranged to perform any of the various methods described herein.

According to yet another aspect of the invention a system is providedthat comprises a plurality of nodes, each node being arranged inaccordance with the above-mentioned apparatus. Each node of theplurality of nodes may be coupled to at least one other node of theplurality of nodes.

According to a further aspect of the invention a computer readablemedium is provided comprising computer readable code operable, in use,to instruct a computer to perform any of the various methods describedherein.

According to another aspect of the invention a method for estimation offrequency of an electrical signal at a node in an electrical systemhaving a plurality of nodes is provided. The method comprisesdetermining a current function indicative of the one or more propertiesof the electrical signal at the node at the current instant in time inaccordance with at least one determined property of the electricalsignal at the node at the current instant in time, a previous functionindicative of the one or more properties of the electrical signal at thenode at a previous instant in time, and a function of the error in theprevious frequency estimation. A current function indicative of one ormore properties of an electrical signal at the at least one other nodeat the current instant in time is received from at least one other nodeof the plurality of nodes. The current function indicative of one ormore properties of the electrical signal at the node and the currentfunction indicative of one or more properties of the electrical signalat the at least one other node are then combined to produce a currentcombined function indicative of one or more properties of the electricalsignal at the node at the current instant in time. Finally, a currentfrequency of the electrical signal at the node is estimated inaccordance with the current combined function.

Also disclosed is a method for estimating a frequency in an electricalpower system having a plurality of nodes, each node being able tomeasure at least one electrical characteristic associated with therespective node and share information related to the at least onemeasured electrical characteristic with at least one other node of theplurality of nodes. The method may comprise measuring, at a node of theplurality of nodes, at least one electrical characteristic associatedwith the node at a current instant in time, k. The method may alsocomprise determining, at the node, a previous function, h(k−m), m=1, 2,. . . , M, indicative of a previous frequency estimation at the node ata previous instant in time, k−m, where typically m=1. Furthermore, themethod comprises determining, at the node, a second function, g(k−m),indicative of a previous frequency estimation at the node at theprevious instant in time, k−m (typically m=1). The second functionincludes a conjugate of the property of the at least one electricalcharacteristic. In addition, the method may comprise determining, at thenode, a function of an error in a previous frequency estimation at thenode. The method may also comprise determining a current function, h,indicative of an estimated frequency of the at least one electricalcharacteristic associated with the node in accordance with the previousfunction, h(k−m), indicative of a previous frequency estimate, thefunction of an error in a previous frequency estimation, and a complexconjugate of the at least one electrical characteristic associated withthe node at the current instant in time, k. Furthermore, the method maycomprise determining a second function, g, of the estimated frequency ofthe at least one electrical characteristic associated with the node inaccordance with the previous second function, g(k−m), indicative of aprevious frequency estimation, the function of an error in the previousfrequency estimation, and the at least one electrical characteristicassociated with the node at the current instant in time, k. In addition,the method may comprise receiving, from at least one other node of theplurality of nodes, the current function indicative of an estimatedfrequency and the current second function of the estimated frequencyassociated with the at least one other node of the plurality of nodes.The method may also comprise utilising, at the node, the currentfunction, h, indicative of an estimated frequency at the node and thecurrent function, h, indicative of an estimated frequency received fromthe at least one other node to produce a current combined (synergetic ordiffused) function indicative of an estimated frequency. Also, themethod may comprise summing, at the node, the current second function,g, indicative of an estimated frequency at the node and the currentfunction, g, indicative of an estimated frequency received from the atleast one other node to produce a current combined (synergetic ordiffused) second function, g, indicative of an estimated frequency.Furthermore, the method may comprise estimating, at each node, a currentfrequency of the at least one electrical characteristic in accordancewith the current combined function and the current combined complexconjugate function.

Embodiments of the invention provide a fast and accurate means forestimating frequency in a three-phase power system.

Embodiments of the invention utilise information collected from aplurality of nodes in a multinode system, which improves the accuracy ofthe frequency estimation.

Embodiments of the invention provide a system having multiplemeasurement nodes, each collecting information indicative of a frequencyof electrical characteristics at each respective node. Each node thenutilises the measurements from neighbouring nodes in order to improvethe frequency estimation accuracy. Each node may weight the informationreceived from each node depending on the perceived quality of themeasurement at that node in order to help improve the accuracy of thefrequency estimation.

Embodiments of the invention relate to a multinode frequency estimationsystem that utilises a widely linear filter in the estimation of thefrequency of an electrical characteristic at the node. The widely linearfilter considers both a function of the electrical measurementcontaining the frequency information and the complex conjugate of thismeasurement.

Embodiments of the invention relate to frequency estimation based onstate space modelling techniques, including Kalman and particlefiltering methods.

A widely linear model may be applied to a Euler representation of theline voltage. In particular, this model may be applied according to:V _(k) e ^(j(ωkT+φ)) =V _(k) cos(ωkT+φ)+jV _(k) sin(ωkT+φ)for which the recursive estimation form is:V _(k) cos(ωkT+φ)=½(V _(k-1) e ^(j(ω(k-1)T+φ)) e ^(jωT) +V _(k-1) e^(−j(ω(k-1)T+φ)) e ^(−jωT))and the corresponding state equation is:

$\begin{bmatrix}x_{k} \\h_{a,k} \\h_{b,k} \\h_{c,k} \\x_{k}^{*} \\h_{a,k}^{*} \\h_{b,k}^{*} \\h_{c,k}^{*}\end{bmatrix} = {\begin{bmatrix}x_{k - 1} \\{h_{a,{k - 1}}x_{k - 1}} \\{h_{b,{k - 1}}x_{k - 1}} \\{h_{c,{k - 1}}x_{k - 1}} \\x_{k - 1}^{*} \\{h_{a,{k - 1}}^{*}x_{k - 1}^{*}} \\{h_{b,{k - 1}}^{*}x_{k - 1}^{*}} \\{h_{c,{k - 1}}^{*}x_{k - 1}^{*}}\end{bmatrix} + u_{k - 1}}$and the observation equation is:

$\begin{bmatrix}v_{a,k} \\v_{b,k} \\v_{c,k}\end{bmatrix} = {{\begin{bmatrix}0 & 1 & 0 & 0 & 0 & 1 & 0 & 0 \\0 & 0 & 1 & 0 & 0 & 0 & 1 & 0 \\0 & 0 & 0 & 1 & 0 & 0 & 0 & 1\end{bmatrix}\begin{bmatrix}x_{k} \\h_{a,k} \\h_{b,k} \\h_{c,k} \\x_{k}^{*} \\h_{a,k}^{*} \\h_{b,k}^{*} \\h_{c,k}^{*}\end{bmatrix}} + n_{k}}$

Where x_(k)=e^(jωT), and for a balanced system:h _(a,k)=½V _(a,k) e ^(j(ωkT+φ))h _(b,k)=½V _(b,k) e ^(j(ωkT+φ−2π/3))h _(c,k)=½V _(c,k) e ^(j(ωkT+φ+2π/3))

The system frequency is the computed as:

${\hat{f}}_{k} = {\frac{1}{2\pi\; T}\;{\arcsin\left( {{??}\left( x_{k} \right)} \right)}}$

A widely linear state space model may be derived using the Euler form sothat the processing is purely on the three-phase voltages rather thanprocessing the Clarke transformed voltages. In addition, the Eulersystem provides improved performance for impulse disturbances, such aslightning. The Euler system may be referred to as a direct widely linearstate space model.

Any dynamic rotating transform may firstly be applied to the Clarke'stransform system voltage. Furthermore, any dynamic rotating transformmay be applied directly to the three system voltages. The dynamicrotating transform may be one the DQ transform or the Park transform.The DQ transform may be applied directly to the three phase voltages,without going through Clarke's transform.

The quantity to be estimated may be the phase Φ of a three-phaseunbalanced system. The phase Φ may be added to the “frequency term” ωkTso that the argument is always (ωkT+Φ)). It may be assumed that thefrequency is changing and the phase is not. However, the phase maychange too, but perhaps not at the same scale. A case may therefore bemade for phase estimation using the same methodology. Such an approachmay be applicable for use with the DQ transformation, where thefrequency term ωkT vanishes, as the system rotates at this rotationalspeed, and only the phase angle Φ is changing.

BRIEF DESCRIPTION OF THE DRAWINGS

Exemplary embodiments of the invention shall now be described withreference to the drawings in which:

FIG. 1 is a simplified diagram of the transmissions and distributionchannels of a power grid system 100;

FIG. 2 is a simplified diagram of a integrated frequency determinationnetwork; and

FIG. 3 shows a structure of a node.

Throughout the description and the drawings, like reference numeralsrefer to like parts.

SPECIFIC DESCRIPTION

FIG. 1 is a simplified diagram of the transmission and distributionchannels of a power grid system 100.

Power is transmitted from a main power station (not shown) at asubtransmission voltage via a high voltage power line 101 and thesubtransmission voltage is then received at a substation 102. Thesubstation 102 steps the voltage down from the voltage of the highvoltage power line 101 used for transmitting the power generated at thepower station long distances to the substation and other substations, toa local distribution voltage delivered to an area local to thesubstation via a distribution voltage line 103. Dividing line 104therefore marks a boundary between the components of the high-voltagetransmission sector 105 of the system 100 and the local distributionsector 106 of the power grid system 100.

The distribution voltage line 103 splits into multiple lines deliveringpower to different parts of the local area. In FIG. 1, the distributionvoltage line splits into three sub-lines delivering power to threeseparate local grids 107, 108,109. The power delivered by thedistribution voltage line 103 is unevenly split between the three localgrids 107, 108, 109 by 20%, 50% and 30% respectively due to the demandof each local grid. For reasons of simplicity only one local grid 109 isdepicted in detail. However, all of the local grids provide one or moreloads that take power from the system and one or more generators thatput power back into the system. Each load may be a particular building,such as a house, apartment block, office or factory, and each generatormay be a local wind farm, household solar PV cell or any other suitablesmall-scale energy harvesting mechanism.

Local grid 109 comprises a local pole transformer 110 which provides afinal main power voltage step down to the standard mains power to bedelivered to the various end points, which in the UK is 240V AC at 50Hz. The local grid 109 provides three local supply lines 111, 112, 113,which each connect to various loads and/or generators. In this caselocal supply line 111 supplies a load 114 and a generator 115, localsupply line 112 supplies three loads 116, 117, 118 and local supply line113 supplies local supply line 119, 120. Each load L will demanddifferent amounts of power and each generator will deliver differentamounts of power and the delivery of such power will be time varying dueto fluctuations in power generation due to, for example in the case ofsolar PV cells, the variation in sunlight intensity, or for windturbines, the variation in wind velocity.

Since the grid system of FIG. 1 has various inputs of power (i.e.generators) and outputs of power (i.e. loads) distributed throughout thesystem, the supply-demand power balance throughout the system willconstantly be varying. The electrical signal at each node (i.e.generator/load) can be considered a separate electrical signal dependenton electrical signals at nearby nodes, or can be considered a singleelectrical signal that has differing characteristics at different nodes.

Due to the variations in the power balance throughout the system of FIG.1, the system runs the risk of various problems occurring such asimbalances between the generation and load, single and dual phasefaults, dual character of load-supply, harmonics, and transientstability, as previously discussed. However, the grid system 100 of FIG.1 is specifically designed to minimise the system failures that occurwithin the grid system 100, as discussed below.

Each node is arranged to measure power system parameters such asvoltages and/or currents and to estimate the frequency of the AC voltageat the node in accordance with those measurements. This frequencyestimation can then be used to detect faults, for example by means ofcomplex circularity analysis of the αβ voltage.

The system is then further designed to improve the speed and accuracy offrequency estimation by enabling each node to share informationregarding its own frequency estimation with some or all of theneighbouring nodes, or remote nodes via a communication link. Each nodethen utilises the information received from neighbouring nodes in itsown frequency estimation, as explained in detail below with reference toFIG. 2.

FIG. 2 is a simplified diagram of an integrated frequency estimationnetwork. FIG. 2 relates to specifically to local grid 109 of FIG. 1.However, it will be appreciated that the network could also includeother local grids, but for ease of explanation only local grid 109 isconsidered in this document. The integrated frequency determinationnetwork of FIG. 2 comprises a plurality of interconnected nodes, whereineach node is any position on the local grid 109 of FIG. 1 where ameasurement of the electrical characteristics at that position is taken.In general, each node will correspond to the position of a load or powersource within the grid, but this is not always necessary.

In general, the system of FIG. 2 works as follows. Each node is arrangedto perform a frequency estimation method to determine the instantaneousfrequency of the electrical signal at that node. However, in order toimprove the performance and accuracy of the frequency estimation, eachnode shares information relating to its frequency determination withneighbouring nodes. Each node then uses the information collecteddirectly from its own node along with information collected fromneighbouring nodes regarding previous frequency estimations in order toassist in the next frequency estimation at that node. Hence, theintegrated frequency determination network is a multinodal networkutilising a multi-input frequency determination scheme, based on thepast measurements (e.g. regression or state space approach). Suchinformation sharing improves the accuracy and speed at which thefrequency is determined and robustness against noise, harmonics, andfaults in communication links.

In FIG. 2, the nodes are communicatively coupled via communicationslinks provided via the power lines, or otherwise. For example, systemssuch as broadband over power lines (BPL) may be utilised. However, itwill be appreciated that it is envisaged that in alternative embodimentsof the invention nodes within such a network may be communicativelyconnected in ways other than by using BPL. For example, all nodes withina network may be wirelessly connected with one another so that all nodesshare information with one another. Alternatively, all nodes may beconnected to each other via a direct wired network or via a centralserver, which routes the information between all nodes. Hence, aninternet-based connection system may be utilised.

Each node in the network receives information from different nodes. Forexample, in FIG. 2, node 114 receives information from nodes 119, 116and 115, node 115 only receives information from node 114, node 119receives information from node 114, 116, and 120, node 116 receivesinformation from node 119, 114 and 117, node 117 receives informationfrom node 116 and node 118, while node 118 receives information fromnode 117. Some nodes therefore receive information from more nodes thanothers. The nodes that receive more information are those nodes that areconnected to a greater number of different nodes. However, the receiptof larger amounts of information should assist such nodes in identifyingsystem errors more quickly. Furthermore, the local and distributednature of the system means that each interconnected node is in factpassing information on to nodes that the node itself is not directlyconnected to. Hence, the network and in particular the passing ofinformation through the network can be complex.

The frequency determination performed at each node shall now bedescribed in detail.

At each node in the network a three phase voltage measurement is takenat a particular instance in time, k. Hence, at each node, threetime-domain voltage measurements are obtained v_(a)(k), v_(b)(k) andv_(c)(k). Each of these phase voltages can be represented mathematicallyin discrete time form as shown by equation 1:

$\begin{matrix}{{{v_{a}(k)} = {{V_{a}(k)}{\cos\left( {{\omega\; k\;\Delta\; T} + \Phi} \right)}}}{{v_{b}(k)} = {{V_{b}(k)}{\cos\left( {{\omega\; k\;\Delta\; T} + \Phi - \frac{2\;\pi}{3}} \right)}}}{{v_{c}(k)} = {{V_{c}(k)}{\cos\left( {{\omega\; k\;\Delta\; T} + \Phi + \frac{2\;\pi}{3}} \right)}}}} & {{Equation}\mspace{14mu} 1}\end{matrix}$where V_(a)(k), V_(b)(k), and V_(c)(k) are the voltage amplitudes of thefirst, second and third phase voltage signals respectively, ω is theangular frequency, k is the instant in time of the voltage measurement,ΔT is the sampling period, and Φ is the phase.

It can be seen from equation 1 that each of the three voltages of thethree phase voltage signal has an associated system frequency, and theirphasors rotate in the complex plane according to that frequency.However, as discussed previously, determining the system frequency fromany one of the three voltages does not provide an accurate and robustmeasure of the overall system frequency at that node at the instance intime, k. Instead, it is necessary to determine the frequenciesassociated with all three voltages.

In order to simplify the process of determining the overall frequency ofa voltage at a node at an instance in time, k, Clarke's transform isutilised, which employs the orthogonal αβ0 transformation to map thetime-dependent three-phase voltage into a zero-sequence v₀ and thedirect- and quadrature-axis components v_(α) and v_(β).

$\begin{matrix}{\begin{bmatrix}{v_{0}(k)} \\{v_{\alpha}(k)} \\{v_{\beta}(k)}\end{bmatrix} = {{\sqrt{\frac{2}{3}}\begin{bmatrix}\frac{\sqrt{2}}{2} & \frac{\sqrt{2}}{2} & \frac{\sqrt{2}}{2} \\1 & {- \frac{1}{2}} & {- \frac{1}{2}} \\0 & \frac{\sqrt{3}}{2} & {- \frac{\sqrt{3}}{2}}\end{bmatrix}}\begin{bmatrix}{v_{a}(k)} \\{v_{b}(k)} \\{v_{c}(k)}\end{bmatrix}}} & {{Equation}\mspace{14mu} 2}\end{matrix}$

For a balanced system, V_(a)(k)=V_(b)(k)=V_(c)(k) and thus v₀(k)=0,v_(α)(k)=A cos(ωkΔT+Φ), and v_(β)(k)=A cos(ωkΔT+Φ+π/2), where v_(α)(k)and v_(β)(k) are the orthogonal coordinates of a point whose position istime variant at a rate proportional to the system frequency. Inpractice, only v_(α) and v_(β) are used and the resulting complexvoltage signal v(k) is given by:v(k)=v _(α)(k)+jv _(β)(k)   Equation 3where j=√{square root over (−1)}.

There is no loss in information in using this representation, and thisvoltage also serves as the desired signal in adaptive frequencyestimation and can be calculated recuresively from:v(k+1)=Ae ^(j((k+1)ωΔT+Φ)) =v(k)e ^(jωΔT)   Equation 4

Where the instantaneous system frequency f is represented by the phasore^(jωΔT)(f=ω/2π).

In normal operating conditions, samples of v(k) are located on a circlein the complex plane with a constant radius A. For a constant samplingfrequency, the probability density function of v(k) is rotationinvariant since v and ve^(jθ) have the same distribution for any real θ.This in turn means that v(k) is second order circular, and in this case,the frequency estimation can be performed adequately by a standardlinear adaptive filter, based on the strictly linear model, on whichtechniques such as the Complex Least Mean Square (CLMS) algorithm arebuilt, and can be summarised as:{circumflex over (v)}(k+1)=v(k)w(k)e(k)=v(k+1)−{circumflex over (v)}(k+1)w(k+1)=w(k)+μp(e(k))q(v*(k))   Equation 5

Where w(k) is the weight coefficient at time instant k, {circumflex over(v)}(k+1) is the estimate of the desired signal v(k+1), e(k) is theestimation error, and μ is the step-size. For the CLMS, the functionsabove take the form p(e(k))=e(k) and q(v(k))=v(k).

Comparing Equation 5 and the strictly linear model in equation 4, thesystem frequency can be estimated from:

$\begin{matrix}{{\hat{f}(k)} = {\frac{1}{2\;\pi\;\Delta\; T}{\sin^{- 1}\left( {{Im}\left( {w(k)} \right)} \right)}}} & {{Equation}\mspace{14mu} 6}\end{matrix}$where Im(⋅) is the imaginary part operator.

While frequency can be estimated from equation 6, such an equation doesnot produce optimal results because it only takes into consideration thesecond order circularity, and is thus only suitable for balancedsystems, for which their Clarke's voltage would be depicted by a circlein a circularity diagram.

As such, in alternative embodiments of the invention alternativefrequency estimation methods are utilised, which produce optimal or nearoptimal results, as discussed below.

When the three-phase power system deviates from its normal condition,such as when the three channel voltages exhibit different levels of dipsor transients, or the phase difference between two channels deviatesfrom the nominal

$\frac{2\;\pi}{3},$v_(a)(k), v_(b)(k) and v_(c)(k) are not identical, and samples of v(k)are not allocated on a circle with a constant radius, the distributionof v(k) is rotation dependent (noncircular), and the signal isaccurately expressed only by using a widely linear model, that is amodel utilising both second-order circular and non-circular information(covariance and pseudocovariance), as defined by:

$\begin{matrix}{\mspace{79mu}{{{v(k)} = {{{A(k)}e^{j{({{\omega\; k\;\Delta\; T} + \phi})}}} + {{B(k)}e^{- {j{({{\omega\; k\;\Delta\; T} + \phi})}}}}}}\mspace{14mu}\mspace{79mu}{Where}\mspace{79mu}{{A(k)} = {\frac{\sqrt{6}\left( {{V_{a}(k)} + {V_{b}(k)} + {V_{c}(k)}} \right)}{6}\mspace{14mu}{and}}}{{B(k)} = {\frac{\sqrt{6}\left( {{2{V_{a}(k)}} - {V_{b}(k)} - {V_{c}(k)}} \right)}{12} - {j\frac{\sqrt{2}\left( {{V_{b}(k)} - {V_{c}(k)}} \right)}{4}}}}}} & {{Equation}\mspace{14mu} 7}\end{matrix}$

In other words, when V_(a)(k), V_(b)(k) and V_(c)(k) are not identical,B(k) is not equal to 0, introducing a certain degree of non-circularity.Since the widely linear model is a second order optimal estimator fornoncircular data, both v(k) and its complex conjugate, v*(k), should beconsidered in the frequency estimation in unbalanced cases. In order toimprove performance, the second order non-circularity (improperness) istherefore also taken into account, i.e.:{circumflex over (v)}(k+1)=v(k)h(k)+v*(k)g(k)   Equation 8

Where h(k) and g(k) are the filter weight coefficients corresponding tothe standard, h, and second, g, updates at time instant k, respectively,wherein the estimation error e(k) is defined as e(k)=v(k+1)−{circumflexover (v)}(k+1).

The weight update of the so-called Augmented Complex Least Mean Square(ACLMS) algorithm, designed for training the widely linear adaptivefilters is given by:h(k+1)=h(k)+μp(e(k))q(v*(k))g(k+1)=g(k)+μq(e(k))q(v(k))   Equation 9where p(e(k))=e(k) and q(v(k))=v(k).

It is noted that there are several stochastic gradient type algorithmsto train the widely linear adaptive filters, such as the normalisedcomplex least mean square (NCLMS), affine projection (AP), least meanphase (LMP), and their sign and leaky versions, depending on the choiceof function used for the functions p(⋅) and q(⋅). The weight update ofall the algorithms use Equation 9 and the specific functions p(⋅) andq(⋅) is summarised as:

$\begin{matrix}{\mspace{79mu}{{{{{NCLMS}\mspace{14mu}{algorithm}\text{:}\mspace{14mu}{p\left( {e(k)} \right)}} = {e(k)}},{{q\left( {v(k)} \right)} = \frac{v(k)}{{{v(k)}}^{2}}}}{{{{AP}\mspace{14mu}{algorithm}\text{:}\mspace{14mu}{p\left( {e(k)} \right)}} = {e(k)}},{{q\left( {V(k)} \right)} = {{V(k)}\left( {{{V^{T}(k)}{V^{*}(k)}} + {\delta\; I}} \right)^{- 1}}}}{{{{LPM}\mspace{14mu}{algorithm}\text{:}\mspace{14mu}{p\left( {e(k)} \right)}} = {< {{v\left( {k + 1} \right)} -} < {\hat{v}\left( {k + 1} \right)}}},{{q\left( {v(k)} \right)} = \frac{j\;{v(k)}}{{\hat{v}}^{*}\left( {k + 1} \right)}}}}} & {{Equation}\mspace{14mu} 10}\end{matrix}$

Where {circumflex over (v)}(k+1) is the at least one property of anelectrical signal determined by estimation, e(k) is the function of anerror in the frequency estimation where defined ase(k)=v(k+1)−{circumflex over (v)}(k+1), V(k) is the input matrix,defined as:

${V(k)} = \begin{bmatrix}{{v(k)},{v\left( {k - 1} \right)},\ldots\mspace{14mu},{v\left( {k - m} \right)}} \\{{v\left( {k - 1} \right)},{v\left( {k - 2} \right)},\ldots\mspace{14mu},{v\left( {k - m - 1} \right)}} \\\vdots \\{{v\left( {k - n} \right)},{v\left( {k - n - 1} \right)},\ldots\mspace{14mu},{v\left( {k - m - n} \right)}}\end{bmatrix}$δ is a regularisation term, which is kept to be a small positiveconstant, I is an identity matrix, ∠ is an angle of a complex variable,and p(⋅) and q(⋅) are functions.

To introduce the corresponding ACLMS-based frequency estimation methodfor the three-phase unbalanced power system, by substituting Equation 7into 8, the estimate {circumflex over (v)}(k+1) can be expressed as:

$\begin{matrix}\begin{matrix}{{\hat{v}\left( {k + 1} \right)} = {{{A(k)}{h(k)}e^{j{({{\omega\; k\;\Delta\; T} + \phi})}}} + {{B(k)}{h(k)}e^{- {j{({{\omega\; k\;\Delta\; T} + \phi})}}}} +}} \\{{{A^{*}(k)}{g(k)}e^{- {j{({{\omega\; k\;\Delta\; T} + \phi})}}}} + {{B^{*}(k)}{g(k)}e^{j{({{\omega\; k\;\Delta\; T} + \phi})}}}} \\{= {{\left( {{{A(k)}{h(k)}} + {{B^{*}(k)}{g(k)}}} \right)e^{j{({{\omega\; k\;\Delta\; T} + \phi})}}} +}} \\{\left( {{{A^{*}(k)}{g(k)}} + {{B(k)}{h(k)}}} \right)e^{- {j{({{\omega\; k\;\Delta\; T} + \phi})}}}}\end{matrix} & {{Equation}\mspace{14mu} 11}\end{matrix}$

While from Equation 7, the expression for v(k+1) can be rewritten as:v(k+1)=A(k+1)e ^(jωΔT) e ^(j(ωkΔT+φ)) +B(k+1)e ^(−jωΔT) e^(−j(ωkΔT+φ))   Equation 12

Therefore, at the steady state, the first term on the right hand side(RHS) of Equation 12 can be estimated approximately by its counterpartin Equation 11; hence, the term e^(jωΔT) containing the frequencyinformation can be estimated as:

$\begin{matrix}{e^{j\;\hat{\omega}\;\Delta\; T} = \frac{{{A(k)}{h(k)}} + {{B^{*}(k)}{g(k)}}}{A\left( {k + 1} \right)}} & {{Equation}\mspace{14mu} 13}\end{matrix}$

Comparing the second term on the RHS of Equation 11 and 12, theevolution of the term e^(−jωΔT) can be expressed as:

$\begin{matrix}{e^{{- j}\;\hat{\omega}\;\Delta\; T} = \frac{{{A^{*}(k)}{g(k)}} + {{B(k)}{h(k)}}}{B\left( {k + 1} \right)}} & {{Equation}\mspace{14mu} 14}\end{matrix}$

Upon taking the complex conjugate, we obtain:

$\begin{matrix}{e^{j\;\hat{\omega}\;\Delta\; T} = \frac{{{A(k)}{g^{*}(k)}} + {{B^{*}(k)}{h^{*}(k)}}}{B^{*}\left( {k + 1} \right)}} & {{Equation}\mspace{14mu} 15}\end{matrix}$

We here adopt the assumptions held implicitly in frequency estimation byadaptive filtering algorithms that A(k+1)≈A(k) and B(k+1)≈B(k), andthus, Equation 13 and Equation 15 can be simplified as:

$\begin{matrix}{e^{j\;\hat{\omega}\;\Delta\; T} = {{h(k)} + {\frac{B^{*}(k)}{A(k)}{g(k)}}}} & {{Equation}\mspace{14mu} 16} \\{e^{j\;\hat{\omega}\;\Delta\; T} = {{h^{*}(k)} + {\frac{A(k)}{B^{*}(k)}{g^{*}(k)}}}} & {{Equation}\mspace{14mu} 17}\end{matrix}$

Since the coefficient A(k) is real-valued, whereas B(k) iscomplex-valued, and thus, B*(k)/A(k)=(B(k)/A(k))*. Since Equation 16should be equal to Equation 17, using a(k)=(B(k)/A(k))*, we can find theform of a(k) by solving the following quadratic equation withcomplex-valued coefficients:g(k)a ²(k)+(h(k)−h*(k))a(k)−g*(k)=0   Equation 18

Where the discriminant of this quadratic equation is given by:

$\begin{matrix}\begin{matrix}{\Delta = \sqrt{\left( {{h(k)} - {h^{*}(k)}} \right)^{2} + {4{{g(k)}}^{2}}}} \\{= {2\sqrt{{- {{Im}^{2}\left( {h(k)} \right)}} + {{g(k)}}^{2}}}}\end{matrix} & {{Equation}\mspace{14mu} 19}\end{matrix}$

Since a(k) is complex-valued, the discriminant is negative, and the tworoots can be found as:

$\begin{matrix}{{{a_{1}(k)} = \frac{{{- j}\;{{??}\left( {h(k)} \right)}} + {j\sqrt{{{Im}^{2}\left( {h(k)} \right)} - {{g(k)}}^{2}}}}{g(k)}}{{a_{2}(k)} = \frac{{{- j}\;{{??}\left( {h(k)} \right)}} - {j\sqrt{{{Im}^{2}\left( {h(k)} \right)} - {{g(k)}}^{2}}}}{g(k)}}} & {{Equation}\mspace{14mu} 20}\end{matrix}$

From Equation 16, the phasor e^(j{circumflex over (ω)}ΔT) is estimatedeither by using h(k)+a₁(k)g(k) or h(k)+a₂(k)g(k). Since the systemfrequency is far smaller than the sampling frequency, the imaginary partof e^(j{circumflex over (ω)}ΔT) is always positive, thus excluding thesolution based on a₂(k). The system frequency {circumflex over (f)}(k)is therefore estimated in the form:

$\begin{matrix}{{\hat{f}(k)} = {\frac{1}{2\;\pi\;\Delta\; T}{\sin^{- 1}\left( {{Im}\left( {{h(k)} + {{a_{1}(k)}{g(k)}}} \right)} \right)}}} & {{Equation}\mspace{14mu} 21}\end{matrix}$

It is noted that, in addition to equation 21, there are severalfunctions to obtain the estimated system frequency, and our aim here isto claim any functions involving h(k−m) and g(k−m), where m can take anyinteger value, that can be used to calculate the frequency estimate.

The above described systems for frequency calculation relate tofrequency calculation at a single node, without interaction between thenodes. The improved global frequency determination shall therefore nowbe discussed in detail.

Whether using the widely linear or state space models, at each node itis possible to derive the functions h and g from the previous error,complex voltage samples, sampling time interval Δ T and the phase Φ. Inparticular, the h and g parameters of each node i are determined fromEquation 9 for the widely linear model or from any one of Equations 22,26, 28, 30 and 32 for the state space models.h _(i)(k+1)= h _(i)(k)+p(e _(i)(k))q(v _(i)*(k))g _(i)(k+1)= g _(i)(k)+p(e _(i)(k))q(v _(i)(k))   Equation 22

Where h_(i)(k) and g_(i)(k) are the h and g estimate at node i at timek, e_(i)(k) is an error calculated using previous error estimates of theparameters h and g together with the signals observed or measured atnode i, namely v_(i)(k−m) where the value of m can be any integer, andfunctions p and q are as in equations 9 and 10. Different methods forestimating h_(i)(k) and g_(i)(k) can be achieved, such as the normalisedleast mean square (NLMS) and least mean phase (LMP), depending on thechoice of function used for the functions p(⋅) and q(⋅). The variables h_(i)(k) and g _(i)(k) are the diffused or global h and g estimates atnode i at time instant k, and are computed as:

$\begin{matrix}{{{\underset{\_}{h}}_{i}(k)} = {{\sum\limits_{l \in N_{i}}{{c\left( {i,l} \right)}{h_{l}(k)}{{\underset{\_}{g}}_{i}(k)}}} = {\sum\limits_{l \in N_{i}}{{c\left( {i,l} \right)}{g_{l}(k)}}}}} & {{Equation}\mspace{14mu} 23}\end{matrix}$

Where all the h estimates received by node i, as well as the h estimatefrom node i itself are weighted using the node dependent weightingfactors or coefficients c(i,l), while the expression N_(i) representsthe neighbourhood of node i, which consists of node i and all the othernodes that can communicate or are connected with node i. It isappreciated that the weighing factors c(i,l) used for computing thediffused variable h _(i)(k) can either be the same as or different fromthose used for g _(i)(k). Moreover, the diffused variables h _(i)(k) andg _(i)(k) can be alternatively computed from the h and g estimatesthrough nonlinear mappings such that:h _(i)(k)=r(h _(N) _(i) )   Equation 24g _(i)(k)=s(g _(N) _(i) )   Equation 25

Wherein r(⋅) and s(⋅) are linear or nonlinear functions, and h_(N) _(i)is a vector consisting of all the h estimates from the neighbourhood ofnode i, similarly, g_(N) _(i) is a vector consisting of all the gestimates from the neighbourhood of node i.

In other words, when each node has determined its local h and g value itthen shares this information with each of the connected nodes in thenetwork. Each node then receives h and g values from any nodes withwhich the respective node is connected. Each node then determines aglobal h and g, based on each of the local h and g values it hasreceived and the h and g value that it has determined.

In practice, each node combines the received h values and received gvalues independently. A set of weighting factors, one associated witheach local node, including itself, are applied to each h and g value inthe summation process. The weighting factors or coefficients c(i,l) areindicative of the perceived accuracy of the information associated witheach node. Hence, some nodes which provide accurate frequency estimationwill be heavily weighted in the summation, while others, which do notprovide accurate frequency estimation, are less heavily weighted in thesummation. A weighting factor of zero could be utilised to discardinformation from a particular node. The node dependent weightings may becontinually updated based upon error determinations.

More specifically, the weighing factors or coefficients c (or r and s)from each node i in the multi-node network can be chosen from a numberof different possible schemes. These schemes include, but are notlimited to:

-   -   Methods for which the weighing factors are chosen according to        the distributed network topology, such as the nearest        neighbourhood method, which involve setting the weighing factors        according to the number of nodes each node can communicate with;    -   Methods involving setting weighing factors according to errors,        be they functions of the estimated or calculated instantaneous        errors, or predetermined errors; or    -   Methods involving using predetermined weighing factors without        any estimation or calculation.

In an alternative embodiment of the invention the variables, h_(i) andg_(i) are estimated using state space modelling techniques, includingKalman and particle filters. Within such frameworks, the estimationproblem is described by state space models, which consists of state andobservation equations, after which the actual algorithm is implementedusing the state space model. The state equation describes how thevariables to be estimated (h_(i) and g_(i)) evolve with time, while theobservation equation describes how the measurements (v_(i)) are relatedto the state variables (h_(i) and g_(i)).

Any one of the following state space models could be utilised for thedetermination of the variables h_(i) and g_(i).

State Space Model 1:

State equation:

$\begin{matrix}{\begin{bmatrix}{h_{i}(k)} \\{g_{i}(k)} \\{h_{i}^{*}(k)} \\{g_{i}^{*}(k)}\end{bmatrix} = {\begin{bmatrix}{h_{i}\left( {k - m} \right)} \\{g_{i}\left( {k - m} \right)} \\{h_{i}^{*}\left( {k - m} \right)} \\{g_{i}^{*}\left( {k - m} \right)}\end{bmatrix} + {u_{i}(k)}}} & {{Equation}\mspace{14mu} 26}\end{matrix}$

Observation equation for a node i in a network:

$\begin{matrix}{\begin{bmatrix}{v_{i}(k)} \\{v_{i}^{*}(k)}\end{bmatrix} = {\begin{bmatrix}{{{v_{i}\left( {k - m} \right)}{h_{i}(k)}} + {{v_{i}^{*}\left( {k - m} \right)}{g_{i}(k)}}} \\{{{v_{i}^{*}\left( {k - m} \right)}{h_{i}^{*}(k)}} + {{v_{i}\left( {k - m} \right)}{g_{i}^{*}(k)}}}\end{bmatrix} + {n_{i}(k)}}} & {{Equation}\mspace{14mu} 27}\end{matrix}$

Where v_(i)(k) is Clarke's αβ transformed voltage at node i at timeinstant k, the variables u_(i)(k) and n_(i)(k) are noise models withpossibly varying statistics, which can be used to model inaccuracies inthe system, as well as to set the convergence rates of algorithms(including Kalman filters), while the symbol (.)* is the complexconjugate operator. The time index variable m can take multiple valuesbut is typically set to 1.

State Space Model 2:

State equation:

$\begin{matrix}{\begin{bmatrix}{h_{i}(k)} \\{g_{i}(k)} \\{z_{i}(k)} \\{h_{i}^{*}(k)} \\{g_{i}^{*}(k)} \\{z_{i}^{*}(k)}\end{bmatrix} = {\begin{bmatrix}{h_{i}\left( {k - m} \right)} \\{g_{i}\left( {k - m} \right)} \\{{{z_{i}\left( {k - m} \right)}{h_{i}\left( {k - m} \right)}} + {{z_{i}^{*}\left( {k - m} \right)}{g_{i}\left( {k - m} \right)}}} \\{h_{i}^{*}\left( {k - m} \right)} \\{g_{i}^{*}\left( {k - m} \right)} \\{{{z_{i}^{*}\left( {k - m} \right)}{h_{i}^{*}\left( {k - m} \right)}} + {{z_{i}\left( {k - m} \right)}{g_{i}^{*}\left( {k - m} \right)}}}\end{bmatrix} + {u_{i}(k)}}} & {{Equation}\mspace{14mu} 28}\end{matrix}$

Observation equation for a node i in a network:

$\begin{matrix}{\begin{bmatrix}{v_{i}(k)} \\{v_{i}^{*}(k)}\end{bmatrix} = {\begin{bmatrix}{z_{i}(k)} \\{z_{i}^{*}(k)}\end{bmatrix} + {n_{i}(k)}}} & {{Equation}\mspace{14mu} 29}\end{matrix}$

The state and observation noise variables, that is variables u_(i)(k)and n_(i)(k), in State Space Model 2 are not necessarily the same asthose in State Space Model 1, due to the differing natures of themodels. State Space Model 2 contains the variable z_(i)(k), which isused to estimate Clarke's transformed voltage v_(i)(k). Consequently,this makes State Space Model 2 more robust to systems with increasedobservation noise.

State Space Model 3:

State equation:

$\begin{matrix}{\begin{bmatrix}{h_{i}(k)} \\{g_{i}(k)}\end{bmatrix} = {\begin{bmatrix}{h_{i}\left( {k - m} \right)} \\{g_{i}\left( {k - m} \right)}\end{bmatrix} + {u_{i}(k)}}} & {{Equation}\mspace{14mu} 30}\end{matrix}$

Observation equation for a node i in a network:v _(i)(k)=v _(i)(k−m)h _(i)(k)+v _(i)*(k−m)g _(i)(k)+n _(i)(k)  Equation 31

State Space Model 3 is essentially the strictly linear version of StateSpace Model 1, and is more suited to cases where the state andobservation noises have circular distributions.

State Space Model 4:

State equation:

$\begin{matrix}{\begin{bmatrix}{h_{i}(k)} \\{g_{i}(k)} \\{z_{i}(k)}\end{bmatrix} = {\begin{bmatrix}{h_{i}\left( {k - m} \right)} \\{g_{i}\left( {k - m} \right)} \\{{{z_{i}\left( {k - m} \right)}{h_{i}\left( {k - m} \right)}} + {{z_{i}^{*}\left( {k - m} \right)}{g_{i}\left( {k - m} \right)}}}\end{bmatrix} + {u_{i}(k)}}} & {{Equation}\mspace{14mu} 32}\end{matrix}$

Observation equation for a node i in a network:v _(i)(k)=z _(i)(k)+n _(i)(k)   Equation 33

State Space Model 4 is the strictly linear version of State Space Model2, and is more suited to systems with circular noise distributions.

The statistics of the noise models, u_(i)(k) and n_(i)(k), in the statespace models above can be different for each state space models, andvary with time. These statistics do not necessarily need to conform tothe true underlying statistics, and can be chosen, computed or estimatedusing a number of different methods. For example, the covariance matrixfor n_(i)(k) can be set to a multiple of the identity matrix (e.g.100I), while simultaneously the covariance matrix for u_(i)(k) is set toa zero matrix or a multiple of the identity matrix (e.g. 0.1I). In fact,the statistics used for these noise models, as well as that for thestate estimate itself, can be varied according to calculated orestimated errors.

The initialisation for all the algorithms based on the state spacemodels and gradient decent techniques described above, are such that theh_(i) and g_(i) variables at each node i are set using the bestestimates available for these variables at each node, which can be basedon the expected nominal operating frequency at each node. However, it isappreciated that the h_(i) and g_(i) variables can be initialised usinga number of different methods, including randomly initialisation, usingpredetermined values or estimated (or computed) prior to starting thealgorithms.

Once the global h and g functions are determined at each node, it ispossible for each node to estimate the frequency of the voltage at thenode. Firstly, the derivation of the frequency calculation utilised byeach node will be explained.

$\begin{matrix}{{f_{i}(k)} = {\frac{1}{2{\pi\Delta}\; T}{\sin^{- 1}\left( {{Im}\left\{ {{{\underset{\_}{h}}_{i}(k)} + {{{\underset{\_}{a}}_{i}(k)}{{\underset{\_}{g}}_{i}(k)}}} \right\}} \right)}}} & {{Equation}\mspace{14mu} 34}\end{matrix}$

Where ΔT is the sampling period, Im{.} is the imaginary part of acomplex number, and

${{\underset{\_}{a}}_{i}(k)} = {\frac{{{- {jIm}}\left\{ {{\underset{\_}{h}}_{i}(k)} \right\}} + {j\sqrt{\left. {{{Im}^{2}\left\{ {{\underset{\_}{h}}_{i}(k)} \right\}} -} \middle| {{\underset{\_}{g}}_{i}(k)} \right|^{2}}}}{{\underset{\_}{g}}_{i}(k)}.}$

A global frequency estimator, consisting of the all the individualestimates by all the nodes in the network, can then be derived as:

$\begin{matrix}{{{f_{global}(k)} = {\frac{1}{2\pi\; T}{\sin^{- 1}\left( {{Im}\left\{ {{{\underset{\_}{h}}_{global}(k)} + {{{\underset{\_}{\alpha}}_{global}(k)}{{\underset{\_}{g}}_{global}(k)}}} \right\}} \right)}}}\mspace{79mu}{{{\underset{\_}{\alpha}}_{global}(k)} = \frac{\begin{matrix}{{{- {jIm}}\left\{ {{\underset{\_}{h}}_{global}(k)} \right\}} +} \\{j\sqrt{\left. {{{Im}^{2}\left\{ {{\underset{\_}{h}}_{global}(k)} \right\}} -} \middle| {{\underset{\_}{g}}_{global}(k)} \right|^{2}}}\end{matrix}}{{\underset{\_}{g}}_{global}(k)}}{{{Where}\mspace{14mu}{{\underset{\_}{h}}_{global}(k)}} = {{\frac{1}{L}{\sum\limits_{i = 1}^{L}\;{{{\underset{\_}{h}}_{i}(k)}\mspace{14mu}{and}\mspace{14mu}{{\underset{\_}{g}}_{global}(k)}}}} = {\frac{1}{L}{\sum\limits_{i = 1}^{K}\;{{\underset{\_}{g}}_{i}(k)}}}}}} & {{Equation}\mspace{14mu} 35}\end{matrix}$

It is noted that, similarly to equation 21, in equations 34 and 35 thereare several functions to obtain the estimated system frequency, and ouraim here is to claim any functions involving h(k−m) and g(k−m), where mcan take any integer value, that can be used to calculate the frequencyestimate.

It will be appreciated that the above-mentioned method for utilisingresources from multiple nodes in order to provide an improved frequencyestimation could equally utilise a strictly linear model, as discussedpreviously. When such a strictly linear model is used the weightcoefficient received from each node is combined to obtain a combinedfunction. The combined function may be determined in a similar way tothe methods described above. Finally, the frequency can then bedetermined based on this combined function.

The structure of the frequency estimation apparatus at example node 200shall now be described with reference to FIG. 3.

Node 200 receives power from an supply power cable 201 at power input202, and the received power is then output to the local area, which maybe a building or such like, via the power output 203 and the local powerdistribution line 204. Sensor 205 measures a voltage at the node 200 andsupplies this node to processor 206, which performs the frequencyestimation. Memory 207 is provided for use during the frequencyestimation process. A communication unit 208 is then coupled to theprocessor and to the internal power connection 209 in order to transmitfrequency estimation information at node 208 to neighbouring nodes, andreceive frequency estimation information from neighbouring nodes.

While the frequency estimation at node 200 has been depicted as astandalone unit, it will be appreciated that such functionality could beincorporated within a standard smart meter. Such smart meters includethe processor, memory, sensor, and communication unit required by thefrequency estimation apparatus. Consequently, the frequency estimationprocess could be incorporated in software downloaded onto a smart meter.The smart meter hardware could then implement the frequency estimationsoftware. Hence, the various methods described above may be implementedby a computer program, as discussed below.

The computer program may include computer code arranged to instruct acomputer to perform the functions of one or more of the various methodsdescribed above. A smart meter may be considered a computer in suchcircumstances. The computer program and/or the code for performing suchmethods may be provided to an apparatus, such as a computer or smartmeter, on a computer readable medium. The computer readable medium couldbe, for example, an electronic, magnetic, optical, electromagnetic,infrared, or semiconductor system, or a propagation medium for datatransmission, for example for downloading the code over the Internet.Non-limiting examples of a physical computer readable medium includesemiconductor or solid state memory, magnetic tape, a removable computerdiskette, a random access memory (RAM), a read-only memory (ROM), arigid magnetic disc, and an optical disk, such as a CD-ROM, CD-R/W orDVD.

An apparatus such as a computer may be configured in accordance withsuch computer code to perform one or more processes in accordance withthe various methods discussed above.

It will be appreciated that any of the state space models or regressionsmodels described above may be utilised for estimating the frequency ofany signal and not just an electrical signal. Furthermore, it will beappreciated that any of the state space models and regression modelsbased on equation (10) described above can be utilised for estimatingfrequency at a single node without receiving and utilising informationfrom neighbouring nodes.

The invention claimed is:
 1. A method for estimating a frequency of anelectrical signal at a node in an electrical system having a pluralityof nodes, the method comprising: determining, at a node of an electricalsystem having a plurality of nodes, at least one property of anelectrical signal at the node at a current instant in time; determininga previous function indicative of one or more properties of theelectrical signal at the node at a previous instant in time; determininga function of an error in a previous frequency estimation at the node ata previous instant in time; determining a current function indicative ofthe one or more properties of the electrical signal at the node at thecurrent instant in time in accordance with the at least one determinedproperty of the electrical signal at the node at the current instant intime, the previous function indicative of the one or more properties ofthe electrical signal at the node at a previous instant in time, and thefunction of the error in the previous frequency estimation; receiving,from at least one other node of the plurality of nodes, a currentfunction indicative of one or more properties of an electrical signal atthe at least one other node at the current instant in time; combiningthe current function indicative of one or more properties of theelectrical signal at the node and the current function indicative of oneor more properties of the electrical signal at the at least one othernode to produce a current combined function indicative of one or moreproperties of the electrical signal at the node at the current instantin time; and estimating a current frequency of the electrical signal atthe node in accordance with the current combined function.
 2. The methodaccording to claim 1, wherein the one or more properties of theelectrical signal that the previous and current functions are indicativeof includes a frequency of the electrical signal.
 3. The methodaccording to claim 2, wherein the one or more properties of theelectrical signal that the previous and current functions are indicativeof includes a phase associated with the electrical signal.
 4. The methodaccording to claim 1, wherein the determined property of the electricalsignal is a voltage, the method further comprising converting thethree-phase electrical signal into a transform-domain, wherein thecurrent function indicative of one or more properties of the electricalsignal in the transform-domain.
 5. The method according to claim 1,wherein the function of an error in the previous estimated frequency isdetermined by multiplication of a calculated error in a previousestimated frequency by a weighting, wherein the weighting is indicativeof a level of confidence in the calculated error.
 6. The methodaccording to claim 1, wherein the function of an error is a function ofa phase-only error.
 7. The method according to claim 1, wherein thecombining is a linear combination process in the original or a transformdomain.
 8. The method according to claim 1, further comprising: afterreceiving the current function indicative of one or more properties ofthe electrical signal at the at least one other node, applying a nodeweighting factor to the received current function before performing thestep of combining, wherein the node weighting factor is indicative of aconfidence in the accuracy of the received current function.
 9. Themethod according to 1, wherein the current estimated frequency{circumflex over (f)}_(i)(k) at the node i, is determined by:${{\hat{f}}_{i}(k)} = {\frac{1}{2{\pi\Delta}\; T}{\sin^{- 1}\left( {{funct}(k)} \right)}}$wherein ΔT is a sampling interval, and funct(k) is the current combinedfunction at the current instant in time, k.
 10. The method according toclaim 1, wherein the previous function indicative of one or moreproperties of the electrical signal at the node at a previous instant intime, the current function indicative of one or more properties of theelectrical signal at the node at the current instant in time, and thefunction received from the at least one other node indicative of one ormore properties of the electrical signal at the at least one other nodeat the current instant in time, each comprise: a first function and asecond function, wherein the first function utilizes the at least onedetermined property of the electrical signal, and the second functionutilizes the complex conjugate of the at least one determined propertyof the signal; and the step of combining comprises: combining, at thenode, the first function associated with the node with the firstfunction associated with the at least one other node to produce acurrent combined first function; and combining, at the node, the secondfunction associated with the node with the second function associatedwith the at least one other node to produce a current combined secondfunction; and the step of estimating comprises estimating a currentfrequency of the at least one electrical characteristic at the node inaccordance with the current combined first function and the currentcombined second function.
 11. The method according to claim 10, whereinthe first function and the second function also utilize a function of acalculated error.
 12. The method according to claim 10, wherein thefirst function and the second function also utilize a function of theconjugate of a complex transform of the at least one determined propertyof the electrical signal at the node at the previous instant in time.13. The method according to claim 10, wherein the first function, h, atthe current instant in time, k, is determined by:h(k)=h(k−m)+μp(e(k−m))q(v*(k−m)) wherein m can take any integer value,h(k−m) is a previous first function, p(e(k−m)) is a function of acalculated error, μ is a weighting factor applied to the calculatederror, and q(v*(k−m)) is a function of the conjugate of a complextransform of the at least one determined property of the electricalsignal at the node at the previous instant in time.
 14. The methodaccording to claim 10, wherein the second function, g, at the currentinstant in time, k, is determined by:g(k)=g(k−m)+μp(e(k−m))q(v(k−m)) wherein m can take any integer value,g(k−m) is a previous second function, p(e(k−m)) is a function of thecalculated error, μ is a weighting factor applied to the calculatederror at the previous instant in time, and q(v(k−m)) is a function of acomplex transform of the at least one determined property of theelectrical signal at the node at a previous instant in time.
 15. Themethod according to claim 10, wherein the current function indicative ofthe one or more properties of the electrical signal at the node at thecurrent instant in time is determined in accordance with the followingstate space model: State equation: $\begin{bmatrix}{h_{i}(k)} \\{g_{i}(k)} \\{h_{i}^{*}(k)} \\{g_{i}^{*}(k)}\end{bmatrix} = {\begin{bmatrix}{h_{i}\left( {k - m} \right)} \\{g_{i}\left( {k - m} \right)} \\{h_{i}^{*}\left( {k - m} \right)} \\{g_{i}^{*}\left( {k - m} \right)}\end{bmatrix} + {u_{i}(k)}}$ Observation equation for the node, i:$\begin{bmatrix}{v_{i}(k)} \\{v_{i}^{*}(k)}\end{bmatrix} = {\begin{bmatrix}{{{v_{i}\left( {k - m} \right)}{h_{i}(k)}} + {{v_{i}^{*}\left( {k - m} \right)}{g_{i}(k)}}} \\{{{v_{i}^{*}\left( {k - m} \right)}{h_{i}^{*}(k)}} + {{v_{i}\left( {k - m} \right)}{g_{i}^{*}(k)}}}\end{bmatrix} + {n_{i}(k)}}$ Wherein h_(i)(k−m) is the first functionindicative of the one or more properties of the electrical signal at thenode, i, at a previous instant in time, k−m, g_(i)(k−m) is the secondfunction indicative of the one or more properties of the electricalsignal at the node at the previous instant in time, k−m, u_(i)(k) andn_(i)(k) are functions indicative of an error in the previous frequencyestimation, and v_(i)(k−m) is the determined property of the electricalsignal at the previous instant in time, k−m, at the node, i.
 16. Themethod according to claim 10, wherein the current function indicative ofthe one or more properties of the electrical signal at the node, i, atthe current instant in time is determined in accordance with thefollowing state space model: $\begin{bmatrix}{h_{i}(k)} \\{g_{i}(k)} \\{z_{i}(k)} \\{h_{i}^{*}(k)} \\{g_{i}^{*}(k)} \\{z_{i}^{*}(k)}\end{bmatrix} = {\begin{bmatrix}{h_{i}\left( {k - m} \right)} \\{g_{i}\left( {k - m} \right)} \\{{{z_{i}\left( {k - m} \right)}{h_{i}\left( {k - m} \right)}} + {{z_{i}^{*}\left( {k - m} \right)}{g_{i}\left( {k - m} \right)}}} \\{h_{i}^{*}\left( {k - m} \right)} \\{g_{i}^{*}\left( {k - m} \right)} \\{{{z_{i}^{*}\left( {k - m} \right)}{h_{i}^{*}\left( {k - m} \right)}} + {{z_{i}\left( {k - m} \right)}{g_{i}^{*}\left( {k - m} \right)}}}\end{bmatrix} + {u_{i}(k)}}$ Observation equation for a node i:$\begin{bmatrix}{v_{i}(k)} \\{v_{i}^{*}(k)}\end{bmatrix} = {\begin{bmatrix}{z_{i}(k)} \\{z_{i}^{*}(k)}\end{bmatrix} + {n_{i}(k)}}$ Wherein h_(i)(k−m) is the first functionindicative of the one or more properties of the electrical signal at thenode i at a previous instant in time, k−m, g_(i)(k−m) is the secondfunction indicative of the one or more properties of the electricalsignal at the node at the previous instant in time, z_(i)(k−m) is anestimation of the one or more properties of the electrical signal at theprevious instant in time, k−m, u_(i)(k) and n_(i)(k) are functionsindicative of an error in the previous frequency estimation, andv_(i)(k−m) is the determined property of the electrical signal at theprevious instant in time, k−m, at the node, i.
 17. The method accordingto claim 1, wherein the current function indicative of the one or moreproperties of the electrical signal at the node is also determined inaccordance with a conjugate of the at least one determined property ofthe electrical signal.
 18. An apparatus for estimating a frequency, theapparatus comprising: a processing unit arranged to perform the methodof claim
 1. 19. A system comprising: a plurality of nodes, each nodecomprising a processing unit arranged to perform the method of claim 1,wherein each node of the plurality of nodes is coupled to at least oneother node of the plurality of nodes.
 20. A non-transitory computerreadable medium comprising computer readable code operable, in use, toinstruct a computer to perform the method of claim 1.